! *******************************************************************************

subroutine read_parameters()
  use myfiles
  use globals
  use grid, only : IntLength,iPer
  implicit none
  
  write(logunit,*)' -reading parameters'
  call flush(logunit)
  
  open(unit=10,file=trim(adjustl(paramfile)),form='formatted',status='old')
  read(10,*)ICDir
  read(10,*)ICBase
  read(10,*)OutDir
  close(10)

end subroutine read_parameters


! ******************************************************************************* 

subroutine read_ic_header_binary()
  use globals
  use myfiles
  use particles
  implicit none
  integer*4 :: file_handle,error,iSnap
  integer*4 :: NumGas,NumDM,NumStars
  character(len=10) :: CharSnap
  character(len=256) :: SnapNoExt
  logical :: imulti

  ! Look for snapshot as single file; if not there look for the .0 etc extension
  SnapNoExt = trim(ICDir)//trim(adjustl(ICBase))

  ! Test for presence of file - if not there, add on .0 to the filename
  Snapshot = trim(SnapNoExt)

  open(10,file=trim(adjustl(Snapshot)),form='unformatted',status='old',iostat=error)
  if(error > 0) then

     ! Can't open, try .0 instead and try again
     Snapshot = trim(SnapNoExt)//'.0'
     open(10,file=trim(adjustl(Snapshot)),form='unformatted',status='old',iostat=error)
     if(error > 0) then
        call abort('  -cannot find snapshot in single- or multi-file format')
     else
        write(logunit,*)'  -opening snapshot in multi-file format'
        imulti = .TRUE.
        close(10)
     endif
 else
     write(logunit,*)'  -opening snapshot in single-file format'
     imulti = .FALSE.
     close(10)
  endif

  ! Read the header
  open(10,file=trim(Snapshot),form='unformatted',status='old')
  read(10) &
       & header%NumPart_ThisFile,  header%MassTable, header%ExpansionFactor,                   &
       & header%Redshift, header%Flag_Sfr, header%Flag_Feedback, header%NumPart_Total,  &
       & header%Flag_Cooling, header%NumFilesPerSnapshot, header%Boxsize,                        &
       & header%Omega0, header%OmegaLambda, header%HubbleParam,                                    &
       & header%Flag_StellarAge, header%Flag_Metals                                                                  
  close(10)

  write(logunit,*)' uuuu',header%HubbleParam
  call flush(logunit)

  NumIn0   = header%NumPart_ThisFile(1)
  NumIn1   = header%NumPart_ThisFile(2)
  NumIn2   = header%NumPart_ThisFile(3)
  NumIn3   = header%NumPart_ThisFile(4)
  NumIn4   = header%NumPart_ThisFile(5)
  NumIn5   = header%NumPart_ThisFile(6)

end subroutine read_ic_header_binary

! *****************************************************************

subroutine read_ic_binary(SnapNoExt)
  use globals
  use myfiles
  use particles
  implicit none
  integer*4 :: error,iSnap,isnapfile,ipart
  integer*4 :: ioffset_by_file(6),ioffset_by_type
  integer*4, allocatable :: temp1d_int(:)
  real*4, allocatable :: temp1d(:),temp2d(:,:)
  character(len=10) :: CharFile,CharFile2
  character(len=256) :: SnapNoExt
     
  integer*4, save :: n_offset(6) = 0

  Snapshot = SnapNoExt
  open(10,file=trim(Snapshot),form='unformatted',status='old')
  read(10) &
       & header%NumPart_ThisFile,  header%MassTable, header%ExpansionFactor,            &
       & header%Redshift, header%Flag_Sfr, header%Flag_Feedback, header%NumPart_Total,  &
       & header%Flag_Cooling, header%NumFilesPerSnapshot, header%Boxsize,               &
       & header%Omega0, header%OmegaLambda, header%HubbleParam,                         &
       & header%Flag_StellarAge, header%Flag_Metals, header%NumPart_Total_HighWord,     &
       & header%Flag_Entropy_Instead_U, header%IC_x0, header%IC_y0, header%IC_z0, header%IC_r0        


  ! Allocate memory for the particles
  call allocate_particles()

  allocate(temp2d(3,sum(header%NumPart_ThisFile)))
  read(10)temp2d
  
  ioffset_by_type = 1
  
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(1)-1
  if(header%NumPart_ThisFile(1) > 0) In0%Coordinates(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(1)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(2)-1
  if(header%NumPart_ThisFile(2) > 0) In1%Coordinates(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(2)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(3)-1
  if(header%NumPart_ThisFile(3) > 0) In2%Coordinates(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(3)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(4)-1
  if(header%NumPart_ThisFile(4) > 0) In3%Coordinates(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(4)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(5)-1
  if(header%NumPart_ThisFile(5) > 0) In4%Coordinates(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(5)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(6)-1
  if(header%NumPart_ThisFile(6) > 0) In5%Coordinates(:,:) = temp2d(:,istart:iend)
  
  deallocate(temp2d)
  
  allocate(temp2d(3,sum(header%NumPart_ThisFile)))
  read(10)temp2d
  
  ioffset_by_type = 1
  
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(1)-1
  if(header%NumPart_ThisFile(1) > 0) In0%Velocity(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(1)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(2)-1
  if(header%NumPart_ThisFile(2) > 0) In1%Velocity(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(2)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(3)-1
  if(header%NumPart_ThisFile(3) > 0) In2%Velocity(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(3)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(4)-1
  if(header%NumPart_ThisFile(4) > 0) In3%Velocity(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(4)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(5)-1
  if(header%NumPart_ThisFile(5) > 0) In4%Velocity(:,:) = temp2d(:,istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(5)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(6)-1
  if(header%NumPart_ThisFile(6) > 0) In5%Velocity(:,:) = temp2d(:,istart:iend)
    
  deallocate(temp2d)

  read(10) ! SKIP the ParticleIDs

  allocate(temp1d(sum(header%NumPart_ThisFile)))
  read(10)temp1d

  ioffset_by_type = 1
  
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(1)-1
  if(header%NumPart_ThisFile(1) > 0) In0%Mass = temp1d(istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(1)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(2)-1
  if(header%NumPart_ThisFile(2) > 0) In1%Mass = temp1d(istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(2)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(3)-1
  if(header%NumPart_ThisFile(3) > 0) In2%Mass = temp1d(istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(3)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(4)-1
  if(header%NumPart_ThisFile(4) > 0) In3%Mass = temp1d(istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(4)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(5)-1
  if(header%NumPart_ThisFile(5) > 0) In4%Mass = temp1d(istart:iend)
  
  ioffset_by_type = ioffset_by_type + header%NumPart_ThisFile(5)
  istart = ioffset_by_type
  iend   = istart + header%NumPart_ThisFile(6)-1
  if(header%NumPart_ThisFile(6) > 0) In5%Mass = temp1d(istart:iend)

  deallocate(temp1d)

  close(10)

  ! Set the ParticleIDs by hand rather than reading in
  do ipart=1,header%NumPart_ThisFile(1)
     In0%ParticleIDs(ipart) = n_offset(1) + ipart
  enddo
  n_offset(1) = n_offset(1) + header%NumPart_ThisFile(1)
  do ipart=1,header%NumPart_ThisFile(2)
     In1%ParticleIDs(ipart) = sum(header%NumPart_Total(1:1)) + n_offset(2) + ipart
  enddo
  n_offset(2) = n_offset(2) + header%NumPart_ThisFile(2)  
  do ipart=1,header%NumPart_ThisFile(3)
     In2%ParticleIDs(ipart) = sum(header%NumPart_Total(1:2)) + n_offset(3) + ipart
  enddo
  n_offset(3) = n_offset(3) + header%NumPart_ThisFile(3)
  do ipart=1,header%NumPart_ThisFile(4)
     In3%ParticleIDs(ipart) = sum(header%NumPart_Total(1:3)) + n_offset(4) + ipart
  enddo
  n_offset(4) = n_offset(4) + header%NumPart_ThisFile(4)
  do ipart=1,header%NumPart_ThisFile(5)
     In4%ParticleIDs(ipart) = sum(header%NumPart_Total(1:4)) + n_offset(5) + ipart
  enddo
  n_offset(5) = n_offset(5) + header%NumPart_ThisFile(5)
  do ipart=1,header%NumPart_ThisFile(6)
     In5%ParticleIDs(ipart) = sum(header%NumPart_Total(1:5)) + n_offset(6) + ipart
  enddo
  n_offset(6) = n_offset(6) + header%NumPart_ThisFile(6)


end subroutine read_ic_binary


subroutine write_new_ics(outfile)
  use hdf5_wrapper
  use myfiles
  use particles
  implicit none
  integer*4 :: file_handle 
  character(len=256) :: outfile

  ! Create the file
  call hdf5_create_file(file_handle,trim(outfile))

  ! Create groups
  call hdf5_create_group(file_handle,'Header')
  call hdf5_create_group(file_handle,'ICParams')
  call hdf5_create_group(file_handle,'PartType1')
  call hdf5_create_group(file_handle,'PartType2')
  call hdf5_create_group(file_handle,'PartType3')

  ! Write the header
  call hdf5_write_attribute(file_handle,'Header/NumPart_ThisFile',header%NumPart_ThisFile)
  call hdf5_write_attribute(file_handle,'Header/NumPart_Total',header%NumPart_Total)
  call hdf5_write_attribute(file_handle,'Header/NumPart_Total_HighWord',header%NumPart_Total_HighWord)
  call hdf5_write_attribute(file_handle,'Header/MassTable',header%MassTable)
  call hdf5_write_attribute(file_handle,'Header/ExpansionFactor',header%ExpansionFactor)
  call hdf5_write_attribute(file_handle,'Header/Redshift',header%Redshift)
  call hdf5_write_attribute(file_handle,'Header/BoxSize',header%BoxSize)
  call hdf5_write_attribute(file_handle,'Header/NumFilesPerSnapshot',header%NumFilesPerSnapshot)
  call hdf5_write_attribute(file_handle,'Header/Omega0',header%Omega0)
  call hdf5_write_attribute(file_handle,'Header/OmegaBaryon',header%OmegaBaryon)
  call hdf5_write_attribute(file_handle,'Header/OmegaLambda',header%OmegaLambda)
  call hdf5_write_attribute(file_handle,'Header/HubbleParam',header%HubbleParam)
  call hdf5_write_attribute(file_handle,'Header/Flag_Sfr',header%Flag_Sfr)
  call hdf5_write_attribute(file_handle,'Header/Flag_Cooling',header%Flag_Cooling)
  call hdf5_write_attribute(file_handle,'Header/Flag_StellarAge',header%Flag_StellarAge)
  call hdf5_write_attribute(file_handle,'Header/Flag_Metals',header%Flag_Metals)
  call hdf5_write_attribute(file_handle,'Header/Flag_Feedback',header%Flag_Feedback)

  ! Write the IC params 
  call hdf5_write_attribute(file_handle,'ICParams/x0',header%IC_x0)
  call hdf5_write_attribute(file_handle,'ICParams/y0',header%IC_y0)
  call hdf5_write_attribute(file_handle,'ICParams/z0',header%IC_z0)
  call hdf5_write_attribute(file_handle,'ICParams/r0',header%IC_r0)

  ! Write the ICs
  call hdf5_write_data(file_handle,'PartType1/Coordinates',In1%Coordinates)
  call hdf5_write_data(file_handle,'PartType2/Coordinates',In2%Coordinates)
  call hdf5_write_data(file_handle,'PartType3/Coordinates',In3%Coordinates)
  call hdf5_write_data(file_handle,'PartType1/Velocity',In1%Velocity)
  call hdf5_write_data(file_handle,'PartType2/Velocity',In2%Velocity)
  call hdf5_write_data(file_handle,'PartType3/Velocity',In3%Velocity)
  call hdf5_write_data(file_handle,'PartType1/ParticleIDs',In1%ParticleIDs)
  call hdf5_write_data(file_handle,'PartType2/ParticleIDs',In2%ParticleIDs)
  call hdf5_write_data(file_handle,'PartType3/ParticleIDs',In3%ParticleIDs)
  call hdf5_write_data(file_handle,'PartType1/Mass',In1%Mass)
  call hdf5_write_data(file_handle,'PartType2/Mass',In2%Mass)
  call hdf5_write_data(file_handle,'PartType3/Mass',In3%Mass)

  ! ******************************************************************************
  ! This section was added to ensure compatibility between these ICs and the HDF5
  ! ICs expected by Marcel van Daalen's customised version of POWMES.
  ! I've just hard wired the OWLS values in here
  ! ******************************************************************************
  ! Units
  call hdf5_create_group(file_handle,'Units')
  call hdf5_write_attribute(file_handle,'Units/UnitDensity_in_cgs',6.76991d-31)
  call hdf5_write_attribute(file_handle,'Units/UnitEnergy_in_cgs',1.989d+53)
  call hdf5_write_attribute(file_handle,'Units/UnitLength_in_cm',3.08568d+24)
  call hdf5_write_attribute(file_handle,'Units/UnitMass_in_g',1.989d+43)
  call hdf5_write_attribute(file_handle,'Units/UnitPressure_in_cgs',6.76991d-21)
  call hdf5_write_attribute(file_handle,'Units/UnitTime_in_s',3.08568d+19)
  call hdf5_write_attribute(file_handle,'Units/UnitVelocity_in_cm_per_s',1.0d5)

  call hdf5_create_group(file_handle,'Constants')
  call hdf5_write_attribute(file_handle,'Constants/AVOGADRO',6.0222d23)
  call hdf5_write_attribute(file_handle,'Constants/BOLTZMANN',1.3806d-16)
  call hdf5_write_attribute(file_handle,'Constants/C',2.9979d10)
  call hdf5_write_attribute(file_handle,'Constants/CM_PER_MPC',3.08568d24)
  call hdf5_write_attribute(file_handle,'Constants/ELECTRONCHARGE',4.8032d-10)
  call hdf5_write_attribute(file_handle,'Constants/ELECTRONMASS',9.10953d-28)
  call hdf5_write_attribute(file_handle,'Constants/EV_TO_ERG',1.60218d-12)
  call hdf5_write_attribute(file_handle,'Constants/GAMMA',1.66667d0)
  call hdf5_write_attribute(file_handle,'Constants/GAS_CONST',8.31425d07)
  call hdf5_write_attribute(file_handle,'Constants/GRAVITY',6.672d-08)
  call hdf5_write_attribute(file_handle,'Constants/HUBBLE',3.24078d-18)
  call hdf5_write_attribute(file_handle,'Constants/PI',3.14159d0)
  call hdf5_write_attribute(file_handle,'Constants/PLANCK',6.6262d-27)
  call hdf5_write_attribute(file_handle,'Constants/PROTONMASS',1.6726d-24)
  call hdf5_write_attribute(file_handle,'Constants/RAD_CONST',7.565d-15)
  call hdf5_write_attribute(file_handle,'Constants/SEC_PER_MEGAYEAR',3.155d13)
  call hdf5_write_attribute(file_handle,'Constants/SEC_PER_YEAR',3.155d07)
  call hdf5_write_attribute(file_handle,'Constants/SOLAR_LUM',3.826d33)
  call hdf5_write_attribute(file_handle,'Constants/SOLAR_MASS',1.989d33)
  call hdf5_write_attribute(file_handle,'Constants/STEFAN',7.5657d-15)
  call hdf5_write_attribute(file_handle,'Constants/THOMPSON',6.65246d-25)
  call hdf5_write_attribute(file_handle,'Constants/T_CMB0',2.728d0)
  call hdf5_write_attribute(file_handle,'Constants/Z_Solar',0.0126637d0)

  call hdf5_close_file(file_handle)

end subroutine write_new_ics
